REPORT  DOCUMENTATION  PAGE 


1.  REPORT  NUMBER 


b.  GOVT  ACCESSI 


AFRL-SR-BL-TR-OO- 


4.  title  (mid  Subti tie)  V  ’  ' 

Advanced  Modeling  of  InP  Crystal  Growth  by  the  MLEK 
Process 


FINAL  —  1  APR  97-31  MAR  98 

6.  PERFORMING  ORG.  REPORT  NUMBER 


7.  AUTHORfiJ 

S.  A.  Orszag,  A.  Tomboulides,  I.  Staroselsky, 
Y.  Zhang,  E.  Barouch 


8.  CONTRACT  OR  GRANT  NUMBERf*; 


F49620-97-C-0011 


9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Cambridge  Hydrodynamics,  Inc. 

P.  O.  Box  1403 
Princeton,  NJ  08542 


II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

AFOSR 

Bolling  AFB,  DC  20332-0001 


10.  PROGRAM  ELEMENT,  PROJECT.  TASK 
area  &  WORK  UNIT  NUMBERS 


12.  report  DATE 


April  1998 


13.  NUMBER  OF  PAGES 


U.  monitoring  agency  name  a  ADDRESS('if  diiierent  from  ControtUng  Office)  15.  SECURITY  CLASS,  fof  reporf; 

Unclassified 

15a.  DECLASSIFICATION/ DOWNGRADING 

schedule 


16.  DISTRIBUTION  STATEMENT  (of  thie  Report) 


Unlimited  -  Unclassified 


17.  DISTRIBUTION  STATEMENT  (of  the  mbatrmct  entered  In  Block  20,  if  different  from  Report) 


Approved  for  public  release, 
distribution  unlimited _ ' 


10.  SUPPLEMENTARY  NOTES 


20000313  031 


19.  KEY  WORDS  (Continue  on  reverae  aide  if  neceaamry  and  identify  by  block  number) 

Crystal  growth  Spectral  methods  Thermal  covection  Rotational  flow  MHD  control 
Boussinesq  equation  Czochralski  technique  Kyropoulos  process 


20.  ABSTRACT  (Continue  on  reverae  aide  If  neceaamry  and  identify  by  block  number) 

Two  and  three-dimensional  flow  of  a  low  Prandtl  number  liquid  metal  is  investigated  including  the 
effects  of  heating,  crucible  and  crystal  rotation,  and  axial  magnetic  fields.  The  crucible  shape  and 
aspect  ratio  are  also  varied  to  determine  their  effects.  A  specially  dsigned  spectral  element  method 
for  axisymmetric  geometries  is  used.  Axisymmetric  flows  that  are  stabilized  with  the  use  of  rotation 
may  be  unstable  to  three-dimensional  disturbances  leading  to  chaotic  behavior.  Axial  magnetic  field 
may  have  a  stabilizing  effect  on  the  flow,  leading  to  axisymmetry;  however,  surprisingly,  large 
magnetic  field  strengths  leads  to  strengthening  of  the  axial  vorticity  and  consequently  strong 
three-dimensionality  and  new  instabilties. 


DD  ,  1473 


imc  QUAUTT INSPSCTBD  3 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Entered) 


Advanced  Modeling  of  InP  Crystal  Growth 
by  the  MLEK  Process 

Final  report  on  AFOSR  Contract  F49620-97-C-001 1  - 


Steven  A.  Orszag 
Ananias  Tomboulides 
Ilya  Staroselsky 
Yansi  Zhang 
Eytan  Barouch 


April  1998 


Cambridge  Hydrodynamics,  Inc. 
P.  O.  Box  1403 
Princeton,  NJ  08542 
(609)  897-9200 


Abstract 


The  two  and  three-dimensional  flow  of  a  low  Prandtl  number  liquid  metal  in  a 
cylindrical  crucible  is  investigated  including  effects  of  heating,  crucible  and  crystal  ro 
tation  and  axial  magnetic  flelds.  The  crucible  shape  and  aspect  ratio  were  varied  to 
study  their  effect  on  the  stability  of  the  flow.  A  spectral  element  numerical  approach, 
specifically  designed  for  axisymmetric  geometries,  is  employed;  this  approach  maintains 
spectral  accuracy  by  removing  the  geometric  singularity  at  the  axis  of  symmetry,  with 
the  use  of  special  Jacobi  polynomials  in  the  radial  direction.  Using  this  approach  in 
the  simulation  of  two-  and  three-dimensional  flows  in  typical  crucible  configurations, 
we  found  that,  as  expected,  reducing  the  crucible  aspect  ratio  and  rounding  its  shape 
results  in  flow  stabilization;  crucible  rotation  and  crystal  differential  rotation  at  speeds 
comparable  to  buoyancy  also  lead  to  flow  stabilization.  However,  axisymmetric  flows 
stabilized  with  the  use  of  rotation,  were  found  to  be  unstable  to  three-dimensional 
disturbances  resulting  in  chaotic  behavior.  Therefore,  solving  the  axisymmetric  prob¬ 
lem  to  find  ranges  of  rotation  speeds  which  stabilize  the  flow  is  not  adequate  for  the 
stabilization  of  3-D  flows.  Moreover,  even  with  strong  crystal  differential  rotation,  the 
resulting  flow  obtained  is  typically  three-dimensional  with  reduced  fluctuation  ampli¬ 
tudes  and  less  chaotic  structure  than  in  the  non-rotating  case.  In  addition,  the  use  of  a 
magnetic  field  aligned  with  the  direction  of  the  crucible  rotation  was  also  investigated 
in  conjunction  with  rotation  and  differential  rotation.  It  was  found  that  the  magnetic 
field  can  have  a  stabihzing  effect  on  the  flow,  and  for  some  strengths  even  renders  the 
flow  axisymmetric.  However,  for  large  values  of  the  magnetic  interaction  parameter  it 
can  lead  to  the  strengthening  of  the  axial  vorticity  component  and,  therefore,  to  more 
three-dimensionahty.  Analysis  of  a  model  problem  involving  Rayleigh-Benard  convec¬ 
tion  in  the  presence  of  rotation  and  magnetic  field  shows  that  the  combined  action  of 
rotation  and  magnetic  field  may  lead  to  the  generation  of  new  instabilities.  More  work 
is  required  in  order  to  understand  the  effect  of  high  values  of  magnetic  field  strength. 

PACS  classification:  81.10,  47.32,  47.65,  02.70.H. 

Keyivords:  Crystal  growth,  spectral  methods,  thermally  driven  rotational  flows, 

MHD  control. 
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1  Introduction 


The  growth  of  silicon,  indium  phosphide  and  other  semi-conductor  crystals  involves  a  number 
of  major  challenges  because  of  the  different  processes  involved,  e.g.  convection  in  the  liquid 
melt  and  the  surrounding  gas,  heat  and  mass  transfer  among  the  system  corhponents  and 
the  sensitivity  of  the  grown  crystals  to  defects,  twinning,  and  other  imperfections.  The 
most  attractive  methods  to  grow  crystals  in  practice  are  the  Czochralski  technique  and  its 
extension  which  involves  use  of  magnetic  fields,  the  magnetic  liquid  encapsulated  Kyropoulos 
(MLEK)  process.  The  complexity  of  these  processes  limits  the  opportunity  to  use  trial-and- 
error  experimental  techniques  to  optimize  crystal  growth.  In  most  cases,  the  flow  in  the  liquid 
melt  is  unsteady  and  three-dimensional  due  to  instabilities  caused  by  thermal  convection, 
and  the  use  of  crucible  and/ or  cystal  rotation,  or  of  magnetic  fields  have  been  shown  to  reduce 
fluctuation  amplitudes.  Thus,  there  is  a  need  for  three-dimensional  numerical  simulations 
which  can  accurately  simulate  the  processes  involved.  Early  3-D  simulations  have  been 
reported  by  Bottaro  k  Zebib  (1989),  Mihelcic  k  Wingerath  (1989),  Jones  (1989),  Leister 
k  Pric  (1992),  Kakimoto  et  al.  (1993),  and  are  mainly  concerned  with  the  onset  of  three- 
dimensionality  and  flow  structure  with  only  thermal  convection.  3-D  melt  flows  during  the 
Czochralski  growth  of  oxide  materials  were  reported  by  Xiao  k  Derby  (1995),  and  eflfects 
of  ampoule  tilting  on  melt  convection  during  Bridgman  growth  were  studied  by  Xiao  et  al. 
(1996).  More  recently,  results  from  3-D  simulations  have  been  reported  at  international 
workshops,  Ben  Hadid  et  al.  (1997),  Tanaka  et  al.  (1997).  Most  numerical  approaches 
are  based  on  finite-difference  or  low-order  finite-element  methods.  Tanaka  et  al.  (1997) 
investigated  pattern  transitions  in  Si  melts  using  a  finite  difference  approach,  for  a  fixed 
value  of  Grashof  number,  Gr,  and  a  range  of  crucible  rotation  speeds.  Ben  Hadid  et  al. 
(1997)  used  a  Legendre  spectral  element  approach  for  the  simulation  of  the  damping  of  3-D 
thermal  convection  of  a  low  Prandtl  number  liquid  in  a  Bridgman  configuration,  for  Gr  of 
the  order  of  10®. 

Here,  we  report  two-  and  three-dimensional  numerical  simulations  of  the  flow  involved 
in  the  MLEK  and  related  growth  processes,  including  numerical  simulations  in  low  Prandtl 
number  melt  convection  with  rotation,  diflFerential  rotation,  and  magnetic  fields.  Our  numer¬ 
ical  approach  is  based  on  spectral  type  methods  (described  below)  specifically  designed  for 
axisymmetric  configurations.  In  the  following  section  we  review  the  conservation  equations 
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and  their  non-dimensionalization  which  gives  rise  to  several  important  non-dimensional  pa¬ 
rameters.  In  section  3  we  go  over  our  numerical  approach;  in  section  4  we  report  results  from 
axisymmetric  simulations  in  different  crucible  configurations,  whereas  in  section  5  we  report 
results  from  full  3-D  simulations,  including  the  effects  of  heating,  rotation  and  differential 
rotation.  Finally,  in  section  6  we  study  the  effect  of  an  axial  magnetic  field  on  the  stability 
of  the  flow  and  make  a  first  attempt  to  explain  our  findings  using  stability  analysis. 


2  Problem  Description 

As  a  first  approximation,  the  crystal-melt  interface  is  assumed  to  be  flat  and  a  typical  con¬ 
figuration  is  shown  in  Fig.  1.  We  have  studied  flows  in  crystal  melts  in  several  situations. 
First,  we  consider  a  model  problem  consisting  of  a  cylindrical  crucible  of  radius  R  and  height 
H  (see  Fig.  1).  Here  we  have  a  prototype  crystal  with  radius  r,  with  r/R  varying  from  0.25 
to  0.5,  being  pulled  from  the  top  surface  of  the  melt.  The  crystal  is  allowed  to  rotate  at  an 
angular  velocity  O2  different  from  that  of  the  crucible  fli.  The  solid  vertical  and  bottom 
boundaries  of  the  crucible  are  maintained  at  a  constant  non-dimensional  temperature  Ti  =  1, 
whereas  the  crystal-melt  interface  is  kept  at  a  fixed  temperature  r2  =  0.  The  top  free  surface 
is  assumed  to  have  negligible  heat  flux  to  the  surrounding  gas  (i.e.  radiation  and  convection 
to  the  gas  is  neglected;  we  note  that  both  assumptions  can  easily  be  removed).  The  equa¬ 
tions  of  motion  are  the  incompressible  Navier-Stokes  equations,  written  in  a  rotating  frame, 
together  with  the  Boussinesq  approximation  for  the  effect  of  the  heating.  The  length  used 
for  non-dimensionalization  is  the  radius  of  the  crucible  R,  whereas  the  non-dimensionalizing 
temperature  is  the  difference  Ti  —  T2. 

The  process  of  crystal  growth  could  potentially  be  represented  by  a  suction  velocity  Vc 
at  the  crystal-melt  interface,  but  this  velocity  is  usually  negligible  when  compared  to  the 
motion  of  the  fluid.  Assuming  that  we  always  work  in  the  frame  of  reference  of  the  rotating 
crucible,  and  that  the  Boussinesq  approximation  is  valid,  the  non-dimensionalized  equations 
of  motion  are 


dv 

-  +  V.VV  = 


— Vp 


k  X  V  ^  1 

2 - I-TH - V^v 

ReE  Re 


(1) 
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(2) 


9T 


1 

RePr 


V^T 


where  k  is  the  direction  of  rotation  (here  aligned  with  the  z  axis),  and  Re  is  equal  to 
Re  =  The  velocity  used  for  the  non-dimensionalization  is  t/  =  {vlR)Gr^l^,  and 

the  non-dimensional  parameters  appearing  in  the  equations  are  the  Grashof  number,  Gr, 
Prandtl  number,  Pr,  and  the  Ekman  number,  defined  as: 


Gr  = 


Pg{T,-T2)R? 


„  u 
Pr  =  — , 
a 


E  = 


u 


R^^li 


(3) 


When  there  is  also  differential  rotation  between  the  crucible  f2i  and  the  crystal  f22,  the 
Rossby  number,  Ro,  is  an  additional  parameter: 


^  2(r2i  -  r^a) 
(ill  -l- 


(4) 


We  note  that  Ro  influences  the  flow  only  through  the  boundary  conditions  at  the  crystal-melt 
interface. 


Typically,  Jones  (1983),  Ristorcelli  &  Lumley  (1992),  heating  of  the  crucible  gives  rise 
to  unsteadiness  and  in  most  cases  three-dimensionality  as  well,  because  natural  convection 
ensues  since  Gr  is  large.  For  most  crystal  growing  applications,  Gr  is  typically  of  0(10^  — 
10*°).  On  the  other  hand,  with  rotation  of  the  crucible,  the  magnitude  of  which  is  described 
by  the  value  of  E,  the  amplitude  of  fluctuations  is  typically  reduced  due  to  centrifugal 
forces  which  are  maximum  along  the  crucible  walls,  which  is  also  where  thermal  convection 
originates.  However,  in  the  presence  of  rotation,  baroclinic  instabilities  can  develop  which 
can  still  render  the  flow  three-dimensional,  Williams  (1971).  To  reduce  heat  flux  fluctuations 
along  the  gas-liquid  interface,  the  crystal  can  be  rotated  independently  of  the  crucible,  with 
strength  determined  by  Ro.  In  this  case,  new  instabilities  can  arise  due  to  the  dynamics  of 
the  shear  layer  between  the  rotating  crystal  and  the  crucible  melt.  In  addition,  rounding  the 
shape  of  the  crucible  to  reduce  sharp  corners,  can  also  reduce  unsteadiness.  The  objective 
is  to  find  the  optimal  rotation  frequencies,  i.e.  E,  Ro,  and  crucible  shapes  for  a  given  Gr 
that  maximally  suppress  thermal  convection  instabilities,  and  lead  to  minimum  fluctuation 
amplitudes  inside  the  melt. 
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3  Numerical  Methodology 


The  numerical  technique  used  in  this  work  is  an  extension  of  the  approach  described  in 
Tomboulides  (1993)  and  Tomboulides  et  al.  (1993).  It  consists  of  a  spectral  element/Fourier 
spatial  discretization  of  the  equations  of  motion  (1),  (2)  in  cylindrical  coordinates,  and  a 
splitting  approach  using  a  second  order  backward  differentiation  scheme  for  time  integration. 
Fourier  decomposition  in  the  azimuthal  direction  gives  the  representation 


M-l 

u{z,  r,  (f>,t)=  Um(r,  z,  (5) 

m=0 

where  m  is  the  azimuthal  wavenumber.  Substituting  (5)  in  the  governing  equations,  and 
applying  the  change  of  variables 


Vm  =  Um  +  iWm,  Wm  =  Vm  -  iWm 


(6) 


gives 


dt 


+  j-^(vvr)  = 


rz  ^2  j  ^ 


dt 


RePr 

dp. 


+  ^„(v.Vv),  =  + 


Re  I,  "  J 


Ur, 


dVm 

dWm 

dt 


■VTm  w  X  V  +  2 


+  X  V  +  2 


k  X 

\  -  - 

(  dPm 

m 

ReE  ) 

1  - 
r 

\  dr 

r 

^  X  vN 

(  dpm 

m 

1 _ 

ReE  ) 

<!> 

[  dr 

n - 

r 

J_  (^1  _  ("^+1)^ 

Re\  r2 


(m  -  1)^ 


Re 


(7a) 
(7b) 
f’m  ("^c) 
U)rri7d) 


where 


+  ii.  (.A) 

dz^  r  dr  \  dr) 


(8) 


^«[  •  lr=^m[  ■  l  +  •  1, 

[  ■  ],^  =  [  ■  ]r  ~  [  ‘  ]</> 
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and  Tin  refers  to  a  Fourier  transform  in  The  coordinate  singularity  at  r  =  0  is  removable, 
since  it  can  be  shown  that  the  behavior  of  the  Fourier  coefficients  of  the  velocity  components 
close  to  the  axis  are 


oc  (/3r”*,7r”‘  (9) 

where  /?  and  7  are  constants,  Orszag  (1974),  Batchelor  (1967).  It  can  be  verified  that 
—  Vm  iwm  is  zero  at  r  =  0  for  all  m  and  scales  like  Vm  +  iwm  oc  a  result 

equivalent  to  the  fact  that  the  vorticity  is  regular  at  r  =  0.  On  the  other  hand,  the  variable 
Wm  =  i^m  —  iWm  has  a  non-zero  value  at  r  =  0  for  m  =  1;  however,  the  coefficient  of  the  1/r^ 
terms  in  (7d)  for  m  =  1  is  zero  and  so  the  singularity  is  removed. 

Numerically,  however,  there  are  still  terms  in  the  equations  where  both  the  numerator  and 
the  denominator  goes  to  zero  at  the  same  rate  close  to  the  axis,  which  means  that  quantities 
of  indeterminate  form  have  to  be  treated.  To  do  this,  a  special  form  of  Jacobi  polynomials 
is  used  as  an  expansion  basis  in  the  r  direction  adjacent  to  the  axis,  which  in  conjunction 
with  L’Hopital’s  rule  results  in  a  removal  of  the  geometrical  singularity,  thus  preserving  the 
spectral  convergence  rate,  Leonard  &  Wray  (1982),  Rpnquist  (1991),  Tomboulides  (1993). 
The  set  of  polynomials  employed  close  to  the  axis  correspond  to  the  Jacobi  polynomials 
p(o.i)^  with  associated  weights  which  are  zero  at  r  =  0. 

For  the  time  integration  of  equations  (7a)-(7d),  we  use  a  fractional  step  method,  in 
conjunction  with  a  mixed  explicit /implicit  stiffly  stable  scheme  of  second  order  of  accuracy 
in  time,  Karniadakis  et  al.  (1991).  A  consistent  Neumann  boundary  condition  is  used  for  the 
pressure,  based  on  the  rotational  form  of  the  viscous  term,  which  nearly  eliminates  splitting 
errors  at  solid  (Dirichlet)  velocity  boundaries,  Tomboulides  et  al.  (1989).  A  total  of  4 
Helmholtz  equations  have  to  be  solved  every  time  step  in  2-D  (5  in  3-D),  one  for  each  of  the 
velocity  components,  and  one  for  the  pressure  and  temperature,  respectively.  The  resulting 
Helmholtz  equations  are  of  the  form 


r  dr  \  dr  j 


Um  ^  '^m  —  Qi 


(10) 


where  Um  stands  for  either  the  velocity,  temperature  or  pressure  azimuthal  Fourier  mode. 
The  constant  is  0  for  the  pressure  equation,  'yoRe/At  for  the  velocity  equations  (70  being 
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a  coefficient  associated  with  the  order  of  the  time-integrating  scheme  used),  and  joR^Pr / 
for  the  temperature  equation. 

The  spatial  discretization  of  Helmholtz  equation  (10)  is  performed  using  two-dimensional 
spectral  elements,  Patera  (1984).  Because  of  the  geometric  singularity  at  r  =  0, 
interpolants  are  used  in  the  radial  direction  in  the  elements  adjacent  to  the  axis  of  symmetry. 
In  elements  not  on  the  axis,  Legendre  interpolants  are  used.  This  kind  of  procedure  is 
necessary  to  achieve  spectrally  accurate  approximations  in  the  whole  computational  domain, 
including  r  =  0.  The  non-linear  terms  are  evaluated  using  collocation  in  physical  space  by 
performing  inverse  Fourier  transforms  on  all  velocity  components  and  temperature,  whereas 
the  rest  of  the  cohiputation  proceeds  in  Fourier  space. 

The  numerical  approach  is  specifically  designed  for  axisymmetric  geometries  and  its  effi¬ 
ciency  relies  on  the  use  of  fast  Fourier  transforms  in  the  azimuthal  direction,  and  fast  banded 
direct  methods  or  conjugate  gradient  based  iterative  methods  in  the  other  two  directions. 
The  resulting  matrix  equations  are  solved  efficiently  by  a  static  condensation  technique  (with 
operation  count  approximately  K^y.  M  x  N^,  where  N  is  the  number  of  grid  points  inside  a 
single  element,  M  is  an  appropriate  bandwidth,  and  Kg  the  total  number  of  elements  (typi¬ 
cally  Ke  <  200  and  <  15).  This  approach  was  used  only  for  the  zeroth  pressure  Fourier 
mode  po  because  of  slow  convergence  properties  when  iterative  techniques  were  used.  The 
Helmholtz  equations  for  the  rest  of  the  unknowns,  i.e.  all  pm,  for  m  ^  0,  and  velocity 
modes,  were  solved  using  preconditioned  conjugate  gradient  iterative  solvers.  The  code  is 
efficiently  parallelized. 


4  Crystal  melt  flow  -  axisymmetric  simulations 

We  report  here  results  from  numerical  experiments  performed  for  crucibles  of  various  aspect 
ratios  and  shapes.  This  section  includes  results  from  2-D  (axisymmetric)  simulations  and 
range  from  flows  with  only  natural  convection,  to  flows  which  include  heating  and  rotation. 
Numerous  resolution  studies  have  been  performed  to  ensure  the  accuracy  and  reliability  of 
the  results  presented  here. 
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4.1  High  aspect  ratio  crucible  H/D=l 

The  Prandtl  number,  Pr,  is  equal  to  0.03,  and  the  Grashof  number,  Gr,  is  kept  equal  to 
2.8  X  10®,  which  is  on  the  low  end  of  the  range  of  Gr  numbers  obtained  in  practice.  Other 
numerical  work  on  buoyancy  driven  flows  in  Si  melts  at  similar  values  of  Gr  was  reported 
by  Mihelcic  k,  Wingerath  (1989).  The  aspect  ratio  of  the  crucible  for  this  set  of  simulations 
is  if/£>  =  1.  The  flow  starts  from  rest  and  Gr  is  gradually  increased  from  10, 000  to  40,  000 
and  then  to  2.8  x  10®.  The  first  two  simulations  reach  a  steady  state  and  then  the  flow 
and  temperature  field  at  Gr  =  40, 000  is  used  as  an  initial  condition  for  the  calculation 
at  Gr  =  2.8  x  10®.  The  time  history  of  several  flow  variables  for  this  run  is  plotted  in 
Fig.  2.  It  is  apparent  that  for  non-dimensional  times  t  <  87.5  even  with  axisymmetry,  the 
flow  is  unsteady  with  moderate  amplitude  fluctuations.  The  result  that,  in  many  flows,  the 
axisymmetric  mode  is  typically  stable,  while  the  first  or  second  azimuthal  Fourier  mode  is 
unstable,  suggests  instability  of  the  three-dimensional  flow  as  well. 

After  the  buoyancy  driven  flow  reaches  a  statistically  steady  state,  rotation  is  introduced 
in  order  to  study  its  effect  on  stabilization  of  the  flow.  The  Ekman  number  is  chosen  to  be 
^  ~  ~  1-67  X  10~^,  which  is  typical  of  rotation  frequencies  in  practice  and  makes 

rotational  effects  of  the  same  order  as  buoyancy  effects.  The  effect  of  rotation  on  the  flow  is 
quite  significant,  since  it  immediately  reduces  the  amplitude  of  the  fluctuations  as  plotted 
in  Fig.  2,  for  times  t  >  87.5. 

The  effect  of  rotation  becomes  more  pronounced  when  isocontours  of  vorticity  are  com¬ 
pared  between  the  two  cases.  The  results  plotted  in  Fig.  3(a)  show  that  without  rotation 
the  flow  develops  two  strong  recirculating  cells  within  which  the  non-dimensional  vorticity 
is  of  order  0(10)  and  the  cold  fluid  is  able  to  convect  all  the  way  to  the  bottom  of  the 
crucible.  On  the  other  hand,  the  results  plotted  in  Fig.  3(b)  show  that  when  rotation  is 
turned  on,  the  two  vortical  cells  disappear  and  the  maximum  non-dimensional  vorticity  in 
the  bulk  of  the  flow  (away  from  walls)  is  reduced  to  order  0(1),  representing  small  amplitude 
fluctuations.  For  both  cases  the  global  maxima  of  vorticity  occur  on  the  side  walls  of  the 
crucible,  as  expected  for  high  Grashof  numbers.  The  fact  that  the  flow  is  stabilized  with 
rotation  in  the  axisymmetric  regime,  does  not  mean  that  the  fluctuation  levels  will  remain 
small  if  the  flow  is  allowed  to  be  three-dimensional  at  these  Gr  numbers.  To  study  this 
question,  three-dimensional  simulations  were  performed  by  using  the  quasi-steady  state  two 
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dimensional  field  with  rotation  at  t  =  172.5  as  an  initial  condition.  These  simulations  are 
described  in  section  5. 

4.2  Effect  of  crucible  aspect  ratio  and  shape 

To  analyze  the  effect  of  crucible  aspect  ratio  on  the  suppression  of  unsteadiness  caused 
by  thermal  convection,  simulations  were  carried  out  for  a  crucible  with  an  aspect  ratio  of 
H/D  =  0.25.  Crucibles  with  aspect  ratios  less  than  1  are  common  in  practice.  The  same 
range  of  Grashof  number  was  investigated  for  the  case  of  a  crucible  with  aspect  ratio  of 
U/D  —  0.25.  Both  2-D  and  3-D  simulations  were  performed  for  this  aspect  ratio,  including 
heating  and  rotation,  however,  only  axisymmetric  simulations  will  be  reported  here.  The  Gr 
number  in  this  case  was  chosen  in  a  way  that  a  comparison  between  the  high  {H/D  =  1)  and 
the  low  {H/D  —  0.25)  aspect  ratio  crucibles  is  meaningful.  Due  to  the  horizontal  as  well  as 
vertical  temperature  gradient,  the  length  scale  used  for  non-dimensionalization  can  be  either 
L  =  {HRy^^  or  L  =  {HR?y^^,  which  results  in  an  equivalent  Gr  number,  for  the  low  aspect 
ratio  crucible,  equal  to  2.8  x  10®  or  4.4  x  10®,  respectively.  This  value  is  equal  or  higher  than 
the  one  used  for  the  high  aspect  ratio  crucible,  and  one  might  expect  a  more  chaotic  flow  to 
be  present.  In  addition,  one  hzis  to  take  into  account  that  now  the  total  volume  of  the  liquid 
metal  contained  in  the  crucible  is  twice  the  one  in  the  high  aspect  ratio  crucible,  with  the 
same  diameter  crystal  grown  on  top  of  it. 

As  can  be  observed  from  the  time  histories  of  flow  components  at  several  points  inside  the 
crucible,  plotted  in  Fig.  4,  the  axisymmetric  flow  in  the  low  aspect  ratio  crucible,  reaches  a 
steady  state  at  long  times.  Comparing  this  result  with  the  2-D  flow  for  the  high  aspect  ratio 
crucible  H/D  =  1,  which  is  strongly  time  dependent  with  0(1)  fluctuations,  it  implies  that 
the  flow  in  the  low  aspect  ratio  crucible  is  more  stable.  Moreover,  isocontours  of  azimuthal 
vorticity  indicate  the  presence  of  only  one  vortical  cell  structure  in  the  crucible,  whereas 
two  such  structures  are  present  in  the  axisymmetric  simulations  for  the  high  aspect  ratio 
crucible.  This  can  be  observed  in  Fig.  5  where  isocontours  of  flow  variables  from  the  steady 
state  field  are  plotted.  As  can  be  observed  in  the  latter  figure,  a  small  pocket  of  stagnant 
fluid,  i.e.  a  small  secondary  recirculation,  exists  close  to  the  lower  right  hand  corner  of  the 
crucible,  which  can  be  a  potential  source  of  instability. 
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As  shown  in  Fig.  7,  this  secondary  recirculation  disappears  when  rounding  the  bottom 
corner  of  the  crucible  in  a  way  that  follows  the  streamlines.  The  flow  again  consists  of  only 
one  recirculating  cell,  but  now  the  crucible  rounded  wall  follows  closely  the  flow  streamline 
structure.  Since  crucibles  used  in  practice  are  rounded  at  the  bottom  without  abrupt  corners, 
our  current  and  future  studies  are  focusing  on  crucible  shapes  similar  to  the  one  in  Fig.  7 
which  is  a  combination  of  a  cylindrical  part  on  top  and  an  ellipsoid  part  at  the  bottom.  The 
flow  for  the  rounded  low-aspect  ratio  crucible  also  reaches  a  steady  state  as  plotted  in  Fig.  6. 


5  Crystal  melt  flow:  3-D  simulations 

In  this  section  we  report  results  from  3-D  simulations  which  range  from  flows  with  only 
natural  convection,  to  flows  which  include  heating,  crucible  rotation  and  crystal  differential 
rotation.  The  long  time  solution  of  the  axisymmetric  flow  for  the  high-aspect  ratio  crucible, 
at  t  =  172.5,  was  used  as  an  initial  condition  for  a  set  of  3-D  simulations.  These  simulations 
started  with  a  perturbation  on  the  first  azimuthal  mode,  rn  =  1,  of  total  energy  of  the  order 
of  10“®.  3-D  simulations  are  reported  here  only  for  the  high  aspect  ratio  crucible,  H/D  =  1. 


5.1  Case  I:  Heating,  Crucible  Rotation,  Crystal  Rotation 

The  growing  crystal  is  usually  rotated  as  it  is  pulled.  The  objective  is  to  improve  uniformity 
by  providing  a  viscous  shear  layer  that  tends  to  isolate  the  growth  interface  from  the  turmoil 
deeper  in  the  melt.  The  crucible  is  also  rotated  to  smooth  out  thermal  asymmetries  that 
might  arise  from  irregularities  in  the  heating.  In  these  simulations,  the  crystal  rotates  in 
the  same  direction  as  the  crucible  but  with  a  Rossby  number  of  Ro  =  2  with  respect  to 
the  crucible.  These  simulations  show  that  the  flow  gradually  develops  three-dimensionality 
which  peaks  in  the  m  —  2  azimuthal  mode,  and  is  caused  by  shear  layer  instabilities  due  to 
the  differential  rotation.  On  the  other  hand,  the  perturbation  on  the  m  =  1  mode  (and  all 
other  odd  modes)  decays  in  time  indicating  absence  of  baroclinic  instability.  In  general,  the 
transition  to  three-dimensionality,  involves  only  even  azimuthal  modes  (rn  =  2,4,6, ..).  The 
time  history  of  the  energy  of  the  first  three  even  azimuthal  modes  are  plotted  as  solid  line 
in  Fig.  8;  the  first  three  odd  modes  are  plotted  as  solid  lines  in  Fig.  9.  A  typical  azimuthal 
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energy  spectrum  is  plotted  in  Fig.  10  where  the  diflFerence  between  the  amplitudes  of  the 
even  versus  odd  modes  is  several  orders  of  magnitude. 

Axial  vorticity  isocontours  are  plotted  in  Figure  (13);  as  can  be  observed  from  this  figure, 
the  flow  has  even  symmetry  with  the  m  =  2  mode  being  the  dominant  one,  and  with  one 
main  cyclonic  vortex  and  two  counter-rotating  ones.  These  results  demonstrate  that  the 
combination  of  rotation  and  differential  rotation  actually  suppress  thermal  convection  fluc¬ 
tuations.  However,  the  dynamics  of  the  shear  layer  created  at  the  crystal-melt  interface  gives 
rise  to  a  new  ordered  flow  structure.  This  structure  is  similar  to  experimental  observations 
reported  in  Hide  &  Titman  (1967).  In  this  set  of  experiments,  involving  isothermal  flow, 
rotation  of  a  disk  inside  a  rotating  crucible  results  in  a  two-  or  higher-fold  flow  symmetry, 
depending  on  the  ratio  of  rotating  speeds  between  the  disk  and  the  crucible.  Therefore,  the 
observed  structure  in  the  melt  flow  is  likely  due  to  the  shear  layer  dynamics  which  dominate 
over  thermal  convection  and  baroclinic  instability. 

Instantaneous  isocontours  of  axial  velocity  and  temperature  are  plotted  in  Figs.  ll(a,b) 
and  12(a,b)  in  the  r—z  and  r—(j)  planes,  respectively,  and  as  can  be  observed  -  especially  from 
the  temperature  isocontours  -  the  flow  demonstrates  an  even  m  =  2  symmetry.  Although 
the  resulting  flow  is  not  axisymmetric,  the  presence  of  more  order  in  the  flow  allows  for  the 
simulation  of  higher  Grashof  numbers  since  the  range  of  excited  scales  is  not  as  large  as  in 
pure  thermal  convection.  For  example,  as  can  be  seen  from  the  results  plotted  in  Fig.  10, 
only  about  16  azimuthal  modes  are  required  to  fully  capture  all  scales  in  the  flow  (to  an 
energy  of  about  10“®). 

5.2  Case  II:  Heating,  Crucible  Rotation 

To  investigate  the  effect  of  crystal  differential  rotation  on  the  stability  of  the  flow,  the  crystal 
rotation  was  turned  off  at  t  =  232.5.  The  time  history  of  azimuthal  mode  energies,  plotted 
with  dashed  lines  in  Figs.  8-9  for  t  >  232.5,  shows  that  three-dimensionality  rapidly  increases 
and,  moreover,  that  all  (both  even  and  odd)  non-axisymmetric  modes  grow  in  time.  The 
amplitude  of  observed  fluctuations  in  this  case  seems  to  be  much  larger  than  for  case  I, 
which  corresponds  to  both  crucible  and  crystal  rotation,  especially  for  the  odd  modes  since 
the  scale  of  the  vertical  axes  in  the  figures  is  logarithmic.  In  addition,  the  flow  structure  for 
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the  case  with  crucible  rotation  only,  plotted  in  Figs.  ll(c,d)  and  12(c,d),  is  more  chaotic 
than  in  case  1,  showing  no  evidence  of  symmetries  in  the  flow. 

5.3  Case  III:  Heating 

The  effect  of  crucible  rotation  was  studied  in  another  numerical  experiment.  This  experiment 
was  performed  by  turning  off  rotation  of  the  crucible  as  well  as  of  the  crystal  at  t  =  80.0, 
and  after  case  I  developed  some  initial  three-dimensionality.  Time  histories  of  azimuthal 
energies  for  this  case,  shown  as  dotted  lines  starting  at  t  =  80.0  in  Figs.  8  and  9,  show  that 
as  soon  as  rotation  is  stopped,  all  non-axisymmetric  modes  increase  sharply  in  energy.  This 
is  an  indication  that  for  the  range  of  parameters  used  in  practice,  pure  thermal  convection 
is  strongly  three-dimensional  and  unsteady,  as  expected,  whereas  rotation  stabilizes  the  flow 
and  reduces  fluctuation  levels.  In  some  cases  it  may  even  be  that  rotation  prevents  three- 
dimensionality.  As  can  be  observed  from  Figs.  ll(e,f)  and  12(e,f)  which  show  isocontours  of 
axial  velocity  and  temperature  on  the  r  —  z  and  r  —  <f>  planes,  the  flow  is  chaotic  and  without 
evidence  of  any  symmetries.  Also,  although  the  Prandtl  number  is  very  low,  low  temperature 
fluid  from  the  top  of  the  crucible  convects  almost  to  the  bottom  of  the  crucible,  indicative 
of  the  large  amplitude  fluctuations  present  in  the  flow.  This  case  was  not  investigated  any 
further  since  the  objective  here  is  not  to  study  the  details  of  turbulent  thermal  convection, 
but  rather  to  identify  ways  to  stabilize  the  flow.  It  should  be  noted  that  simulations  for  case 
III  are  very  costly  in  terms  of  computational  time  and  memory  because  of  the  wide  range  of 
excited  scales. 


6  Magnetic  field  simulations 

In  some  crytal  growing  applications  the  crystal  melt  is  stabilized  by  the  use  of  a  magnetic 
field.  The  use  of  a  magnetic  field,  aligned  with  the  rotation  axis,  for  the  stabilization  of 
the  flow  was  investigated  in  a  series  of  numerical  experiments.  Since  the  magnetic  Reynolds 
number  is  very  low,  the  “induction-less”  approximation  is  used  which  leads  to  a  simplified  set 
of  equations.  The  non-dimensionalized  governing  equations  at  the  low  magnetic  Reynolds 
number  limit  are 
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dv 

— +  v- Vv  =  -Vp- 


2k  XV  ^ 
ReE 


-^VV  +  r  +  iVjxBo 

Re 


j  =  -V^  + V  X  Bo 


VV  =  V  •  (v  X  Bo) 

-=n.(vxB.) 

Here,  Re  is  again  defined  as  Re  —  Bo  is  the  imposed  axial  magnetic  field,  and  (f) 

is  the  electric  potential  corresponding  to  the  induced  electric  field.  The  non-dimensionalized 
magnetic  field  strength.  Bo,  is  equal  to  the  unit  vector  k  in  the  axial  direction.  The  boundary 
conditions  used  for  the  electric  potential  are  vanishing  normal  components  of  the  electric  cur¬ 
rent  density  j  at  the  crucible  walls  as  well  as  at  the  crystal  interface.  The  additional  parame¬ 
ter  which  appears  in  the  equations  is  the  magnetic  interaction  parameter,  N  =  BlRem/pU^- 
N  is  related  to  the  Hartmann  number  by 


Ha'^  =  NRe  =  —^RemRe 
pU^ 

Various  simulations  were  performed  with  the  main  focus  being  the  understanding  of  the 
influence  of  an  axial  magnetic  field  on  the  melt  flow.  All  simulations  were  performed  starting 
from  case  I,  which  involves  crucible  rotation  and  crystal  differential  rotation.  After  case  I 
reaches  a  quasi-steady  state,  the  magnetic  field  is  turned  on  and  its  influence  on  the  flow 
is  monitored  in  time.  In  Fig.  14,  the  time  history  of  azimuthal  mode  energies  is  plotted 
for  the  cases  with  AT  =  0,  1,  and  20.  As  can  be  observed,  at  AT  =  1,  all  non-axisymmetric 
modes  are  damped  resulting  in  an  almost  axisymmetric  flow.  On  the  other  hand,  at  A^  =  20, 
the  amplitude  of  non-axisymmetric  modes  increases.  In  fact,  at  an  intermediate  value,  e.g. 
N  =  5,  the  amplitude  of  the  non-axisymmetric  modes  is  reduced  but  not  to  zero.  In  order  to 
visualize  the  effect  of  the  magnetic  fleld  on  the  flow  field,  streamwise  vorticity  isocontours  at 
z  =  I  are  plotted  in  Figs.  15(a-d).  As  can  be  seen  in  these  figures,  the  resulting  vortex  shape 
for  A  =  1  (Fig.  15(a))  is  almost  a  circle,  corresponding  to  an  axisymmetric  flow;  on  the 
other  hand,  for  A  =  5  (Fig.  15(b)),  the  shape  of  the  main  vortex  is  only  slightly  distorted 
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from  axisymmetry,  whereas  for  =  10  (Fig.  15(c))  and  20  (Fig.  15(d))  the  main  central 
vortex  is  distorted  from  a  circular  shape  to  an  ellipse  with  a  higher  aspect  ratio  than  for 
the  case  without  magnetic  field,  plotted  in  Fig.  13.  In  addition,  the  vortex  structure  rotates 
with  respect  to  the  rotating  crucible,  resulting  in  an  almost  periodic  flow  in  all  cases. 

Our  simulations  indicate  that  the  imposed  axial  magnetic  field,  can  either  dampen  three- 
dimensionality,  or,  in  certain  parameter  regimes,  enhance  it,  and  thereby  destabilize  the 
flow.  In  the  following  subsection,  we  present  a  simplified  model  that  illustrates  this  effect 
of  magnetic  fields  in  the  presence  of  rotation.  In  this  model,  we  have  extended  the  analysis 
in  Chandrasekhar  (1961)  to  analyze  the  effects  of  magnetic  fields  on  the  stability  of  Benard 
convection.  This  analysis  which  involves  the  solution  of  a  12th  order  boundary  value  problem, 
leads  to  interesting  insights  on  the  mechanisms  by  which  magnetic  fields  can  be  destabilizing, 
as  shown  in  the  following  subsection. 

6.1  Influence  of  rotation  and  magnetic  field  on  the  stability  of 
thermal  convection 

Crystal  rotation  and  magnetic  fields  individually  may  suppress  certain  instabilities,  however, 
combined  together  they  may  lead  to  the  generation  of  new  instabilities.  Here,  we  analyze 
an  example  of  such  new  instabilities  which  can  arise  by  the  combined  use  of  rotation  and 
magnetic  field.  We  consider,  following  Chandrasekhar  (1961),  the  case  of  Rayleigh-Benard 
convection  in  the  presence  of  rotation  and  magnetic  field.  As  described  in  section  2,  the  flow 
is  characterized  by  the  Rayleigh  number,  i2a,  Prandtl  number,  Pr,  Taylor,  T,  or  Ekman,  E, 
numbers,  with  the  Taylor  number  being  T  =  1/El^,  and  the  Hartmann  number,  Ha. 

Separately  both  rotation  and  magnetic  field  inhibit  the  onset  of  instability  and  they  both 
elongate  the  cells  which  appear  at  margined  stability.  These  effects  have  a  common  origin; 
the  flow  becomes  more  two-dimensional.  Acting  together,  however,  they  may  lead  to  a 
new  instability.  To  understand  the  origin  of  this  paradoxical  behavior,  we  note  that  rotation 
induces  a  component  of  vorticity  in  the  direction  of  rotation  resulting  in  streamlines  becoming 
closely  wound  spirals  with  motions  principally  confined  to  planes  transverse  to  the  rotation. 
On  the  other  hand,  magnetic  fields  suppress  motions  transverse  to  their  direction  so  motion 
along  magnetic  field  lines  becomes  dominant.  In  addition,  instability  in  the  presence  of 


rotation  alone  is  a  Hopf  bifurcation,  although  it  appears  as  exchange  of  stability  (stationary 
convection)  when  a  magnetic  field  is  present  alone.  Acting  together,  the  conflicting  behavior 
of  the  effects  of  rotation  and  axial  magnetic  field  may  lead  to  a  reinforcement  of  each  other,  i.e. 
instability,  under  certain  conditions.  This  complex  behavior  leading  to  instability  requires 
full  mathematical  analysis;  an  example  is  given  below. 

In  the  simplest  case  of  two  free  boundaries  the  characteristic  equation  for  the  critical 
Rayleigh  number  takes  the  form 

_  ^4  +  Ti(l  +  3;)} 

x[{l+xY  +  Ha?] 

where 

T 

x  =  --,  and  =  -  (12) 

and  £  is  the  characteristic  scale  of  the  onset  of  instability.  These  equations  determine  the 
instability  threshold  in  the  case  of  the  onset  of  instability  as  stationary  convection.  In  the 
case  of  the  onset  of  instability  as  overstability  the  characteristic  equations  take  the  form 


where 


Ra  =  {1+xY  +  Ha^ 
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These  equations  are  accurate  in  the  case  of  low  magnetic  Prandtl  number  which  is  the  case 
for  crystal  growth. 


The  solution  of  these  equations  for  Ekman  number  equal  to  E  =  10“®  (or  Taylor  number 
T  =  10®)  and  Prandtl  number  Pr  =  0.2  is  shown  in  Fig.  (16).  Let  us  consider  the  stability 
of  the  flow  with  increasing  magnetic  field,  starting  from  Rayleigh  numbers  for  which  the  flow 
is  stable  in  the  absence  of  the  magnetic  field.  As  can  be  observed,  the  flow  may  become 
unstable  at  a  certain  level  of  the  magnetic  field.  For  this  model  problem,  if  we  increase  the 
strength  of  the  magnetic  field  even  further  the  flow  stabilizes  once  again. 

The  model  problem  is  meant  to  illustrate  the  complex  combined  effect  of  rotation  and 
magnetic  fields  in  the  case  of  Rayleigh-Benard  convection  in  an  infinite  layer.  In  other  more 
complex  systems,  the  effect  of  the  destabilization  of  thermal  convection  by  the  action  of  a 
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magnetic  field  in  the  presence  of  rotation  is  sensitive  to  the  geometry,  Prandtl  number  and 
configuration  of  temperature  gradients.  More  complex  configurations  cannot  be  analyzed 
using  such  simple  models,  and  the  stability  of  flows  in  crystal  melts  in  the  presence  of 
magnetic  field  and  rotation  can  only  be  studied  using  three-dimensional  direct  numerical 
simulation. 


7  Conclusions 

The  2-D  and  3-D  flow  of  a  low  Prandtl  number  liquid  metal  in  a  cylindrical  crucible  were 
investigated  including  effects  of  heating,  crucible  and  crystal  rotation  and  axial  magnetic 
fields  using  a  spectral  element  numerical  approach,  specifically  designed  for  axisymmetric 
geometries.  It  was  found  that,  as  expected,  reducing  the  crucible  aspect  ratio  and  rounding 
its  shape  results  in  flow  stabilization.  Crucible  rotation  and  crystal  differential  rotation,  at 
speeds  comparable  to  buoyancy,  enhance  flow  stabilization  in  2-D  as  well  as  in  3-D.  However, 
axisymmetric  flows  stabilized  with  the  use  of  rotation,  were  found  to  be  unstable  to  three- 
dimensional  disturbances  resulting  in  chaotic  behavior.  Therefore,  solving  the  axisymmetric 
problem  to  find  ranges  of  rotation  speeds  which  stabilize  the  flow  is  not  adequate  for  the  sta¬ 
bilization  of  3-D  flows.  Moreover,  even  with  strong  crystal  differential  rotation,  the  resulting 
flow  obtained  is  typically  three-dimensional  with  reduced  fluctuation  amplitudes  and  less 
chaotic  structure  compared  to  the  non-rotating  case.  In  addition,  the  use  of  a  magnetic  field 
aligned  with  the  direction  of  the  crucible  rotation  was  also  investigated  in  conjunction  with 
rotation  and  differential  rotation.  It  was  found  that  the  magnetic  field  can  have  a  stabilizing 
effect  on  the  flow,  and  for  some  strengths  even  render  the  flow  axisymmetric.  However,  for 
large  values  of  the  magnetic  interaction  parameter  it  can  lead  to  the  strengthening  of  the 
axial  vorticity  component  and,  therefore,  to  more  three-dimensionality.  Analysis  of  a  model 
problem  involving  Rayleigh- Benard  convection  in  the  presence  of  rotation  and  magnetic  field 
shows  that  the  combined  action  of  rotation  and  magnetic  field  may  lead  to  the  generation 
of  new  instabilities.  More  work  is  required  in  order  to  understand  the  effect  of  high  values 
of  magnetic  field  strength. 
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9  Figure  Captions 

•  Figure  1:  Geometric  configuration  of  a  crystal  melt  crucible. 

•  Figure  2:  Time  history  of  u,v,T  for  Gr  =  2.8  x  10®,  without  rotation  until  t  >  87.5, 

and  with  rotation  with  E  =  =  1.67  x  10“®  after  t  >  87.5.  Axisymmetric  flow 

simulation. 

•  Figure  3:  Instantaneous  vorticity  field  (a)  without  rotation  and  (b)  with  rotation,  for 
a  crucible  with  H/D  =  1.  Axisymmetric  flow  simulation.  Note  that  this  flgure  was 
generated  by  transformation  of  color  plots  to  grey  scale,  so  the  darkest  tone  does  not 
correspond  to  the  highest  value. 

•  Figure  4:  Time  history  of  velocity,  pressure  and  temperature  for  axisymmetric  flow  in 
a  low  aspect  ratio  crucible  at  Gr  =  2.8  x  10®. 

•  Figure  5:  Isocontours  of  (a)  temperature,  and  (b)  vorticity  for  the  steady  state  flow 
field  in  a  crucible  with  aspect  ratio  H/D  —  0.25,  and  Gr  =  2.8  x  10®.  Note  that  this 
figure  was  generated  by  transformation  of  color  plots  to  grey  scale,  so  the  darkest  tone 
does  not  correspond  to  the  highest  value. 

•  Figure  6:  Time  history  of  velocity,  pressure  and  temperature  for  axisymmetric  flow  in 
a  rounded  crucible  at  Gr  =  2.8  x  10®. 

•  Figure  7:  Isocontours  of  (a)  temperature,  and  (b)  vorticity  for  the  steady  state  flow 
field  in  a  rounded  crucible  with  aspect  ratio  H/D  —  0.25,  and  Gr  =  2.8  x  10®.  Note 
that  this  figure  was  generated  by  transformation  of  color  plots  to  grey  scale,  so  the 
darkest  tone  does  not  correspond  to  the  highest  value. 

•  Figure  8:  Time  history  of  even  azimuthal  mode  energies. 

•  Figure  9:  Time  history  of  odd  azimuthal  mode  energies. 

•  Figure  10:  A  typical  azimuthal  energy  spectrum  for  flow  with  rotation  and  differential 
rotation. 

•  Figure  11:  Instantaneous  isocontours  of  axial  velocity  (left)  and  temperature  (right) 
on  a  r  —  z  plane  for  cases  (a),(b)  with  heating  rotation,  and  differential  rotation,  (c), 
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(d)  with  heating  and  rotation  only,  (e),  (f)  with  heating  only.  Note  that  this  figure 
was  generated  by  transformation  of  color  plots  to  grey  scale,  so  the  darkest  tone  does 
not  correspond  to  the  highest  value. 

•  Figure  12:  Instantaneous  isocontours  of  axial  velocity  (left)  and  temperature  (right) 
on  a  r  —  0  plane  for  cases  (a),(b)  with  heating  rotation,  and  differential  rotation,  (c),(d) 
with  heating  and  rotation  only,  (e),(f)  with  heating  only.  All  axial  velocity  isocontours 
are  at  2:  =  1,  whereas  all  temperature  isocontours  are  at  2  =  1.75.  Note  that  this  figure 
was  generated  by  transformation  of  color  plots  to  grey  scale,  so  the  darkest  tone  does 
not  correspond  to  the  highest  value. 

•  Figure  13:  Isocontours  of  axial  vorticity  at  2  =  1  for  fiow  with  heating  rotation,  and 
differential  rotation.  Note  that  this  figure  was  generated  by  transformation  of  color 
plots  to  grey  scale,  so  the  darkest  tone  does  not  correspond  to  the  highest  value. 

•  Figure  14:  Time  history  of  the  even  azimuthal  modal  energies  when  different  strength 
magnetic  fields  axe  applied. 

•  Figure  15:  Isocontours  of  axial  vorticity  at  2  =  1  for  fiow  with  a  magnetic  field  with 
a)  A'^  =  1,  b)  iV  =  5,  c)  A  =  10,  and  d)  N  =  20.  Note  that  this  figure  was  generated 
by  transformation  of  color  plots  to  grey  scale,  so  the  darkest  tone  does  not  correspond 
to  the  highest  value. 

•  Figure  16:  The  stability  diagram  of  critical  Rayleigh  number  as  a  function  of  Hartmann 
number  for  Ekman  number  10“^ 
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